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We consider problems of Bayesian inference for a spatial epidemic on a graph, where the final 
state of the epidemic corresponds to bond percolation, and where only the set or number of finally 
infected sites is observed. We develop appropriate Markov chain Monte Carlo algorithms, demon- 
strating their effectiveness, and we study problems of optimal experimental design. In particular, 
we demonstrate that for lattice-based processes an experiment on a sparsified lattice can yield more 
information on model parameters than one conducted on a complete lattice. We also prove some 
probabilistic results about the behaviour of estimators associated with large infected clusters. 
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I. INTRODUCTION 

Many real-world phenomena can be modelled by ran- 
dom graphs, or more generally, by dynamically chang- 
ing random graphs. Specifically, host-pathogen biolog- 
ical systems that may combine primary and nearest- 
neighbour or long-range secondary infection processes 
can be efficiently described by spatio-temporal models 
based on random graphs evolving in time 1]. 

Although a continuous observation of an epidemic is 
not always possible, a spatial 'snapshot' may provide one 
with some, albeit highly incomplete, knowledge about the 
epidemic. In terms of the model this knowledge results 
in a random graph realised in some metric space. More- 
over, under some circumstances it is not even possible 
to observe some or all of the edges of such a random 
graph — all one would know then are the vertices which 
correspond to the infected sites, i.e. to those sites which 
interacted as a result of the evolution of the process under 
consideration. 

One particular application refers to the colonisation of 
susceptible sites, such as seeds or plants grown on a lat- 
tice, by virus, fungal, or bacterial pathogens with limited 
dispersal abilities. A typical example is the spread of in- 
fections through populations of seedlings by the fungal 
pathogen, Rhizoctonia solani Kiihn. This economically- 
important pathogen is wide spread with a remarkably 
wide host range In addition to its intrinsic economic 
importance, it has been extensively used as an experi- 
mental model system to test epidemiological hypotheses 
in replicated microcosms [5, '3| and to study biological 
control of pathozone behaviour by an antagonistic fungus 
and disease dynamics i5i] . Transmission of infection be- 
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tween plants occurs by mycelial growth from an infected 
host, with preferential spread along soil surfaces — hence 
the missing information about the structure of interac- 
tions. 

The spread of infections with limited dispersal 
abilities among plants can be viewed as a spatial 
SIR (susceptible infective — >■ removed) epidemic with 
nearest-neighbour secondary infections and removals, 
and can be related to a percolation process on a finite 
static network P, 0| . Percolation has also been widely 
used to model various phenomena of disordered media 
in physics, chemistry, biology, engineering and materials 
science 0, II] , and also forest fires 

Bayesian estimation for percolation models of disease 
spread in plant populations in the context of the spread of 
Rhizoctonia solani has been presented in [l[. In the 
spread of this soil-borne fungal plant pathogen among 
discrete sites of nutrient resource is studied using simple 
concepts of percolation theory; a distinction is made be- 
tween invasive and noninvasive saprotrophic spread. The 
authors of |ll, Q formulated statistical methods for fit- 
ting and testing percolation-based, spatio-temporal mod- 
els that are generally applicable to biological or physical 
processes that evolve in time in spatially distributed pop- 
ulations. Estimation of spatial parameters from a single 
snapshot of an epidemic evolving on a discretised grid 
under the assumption that fundamental spatial statistics 
are near equilibrium is studied in [11]. 

The difficulties in performing inference for these mod- 
els in the presence of observational uncertainty or in- 
complete observations can be overcome to an extent by 
employing a Bayesian approach and modern powerful 
computational techniques — mainly Markov chain Monte 
Carlo (McMC), e.g. [lij. McMC methods often offer 
important advantages over existing methods of analysis. 
Particularly, they allow a much greater degree of mod- 
elling flexibility, although the implementation of McMC 
methods can be problematic because of convergence and 
mixing difficulties which arise due to the amount and 
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nature of missing data. 

An aspect which has received httle attention in the 
context of the described models is that of experimental 
design. Statisticians have investigated the question of 
experimental design in the Bayesian framework (see 13 1 
for a review). The work of Miiller and others [3, 15 1 
examined the ways of identifying designs that maximise 
the expectation of a utility function. 

The goal of this paper is to study inference and optimal 
design problems for a model which is obtained by tak- 
ing a final snapshot of nearest-neighbour disease spread 
spatio-temporal dynamics and relating it to the percola- 
tion process on a regular grid. We use the utility-based 
Bayesian framework for this purpose and discuss generic 
issues that arise in this context. 

The paper is arranged as follows. In Section 2 we de- 
scribe the percolation model considered and describe the 
partial sampling scenarios that form the focus of the pa- 
per. In Section 3 we describe how, by embedding the 
problem in a framework of Bayesian inference, it is pos- 
sible to estimate likelihood functions for percolation mod- 
els under the respective sampling scenarios using Markov 
chain Monte carlo (McMC) methods. Moreover we prove 
some theoretical results regarding the limiting properties 
of the likelihood for the case when only cluster size is 
observed, and elicit interesting connections with percola- 
tion thresholds. In Sections 5 and 6 we turn our attention 
to the question of experimental design and, in particu- 
lar, the identification of experimental lattices that can 
maximise the information on the percolation probabil- 
ity yielded by an experimental process observed on it, 
when only the vertex set of the connected cluster is ob- 
served. We demonstrate that the optimal lattice may be 
one which includes regions which have been sparsified by 
removal of certain vertices and edges. These results have 
important implications for the design of experiments on, 
for example, fungal pathogen spread in populations of 
agricultural plants, where hosts are typically arranged 
on a lattice. 



II. NEAREST-NEIGHBOUR INTERACTION 
MODEL AND PERCOLATION 

A. Model description 

Disease spread as a result of (typically) short-range 
contact between, for example, plants can be modelled as 
a transmission process on an undirected graph. Nodes, 
or vertices, of the graph correspond to possible locations 
of plants, and edges of the graph link locations which are 
considered to be neighbours. In a classical SIR model 
each node, or vertex, of the graph is in one of three 
states: either it is occupied by a healthy, but suscepti- 
ble, plant (state 5), or it is occupied by an infected and 
infectious plant (state /), or finally it is empty, any plant 
at that location having died and thus being considered 
removed (state R). A plant at node i, once infected (or 



from time if initially infected), remains in the infected 
state / for some random lifetime after which it dies, 
so that node i then remains in the empty state R ever 
thereafter. During its infected lifetime the plant at node i 
sends further infections to each of its neighbouring nodes 
j as a Poisson process with rate (so that the proba- 
bility that an infection travels from i to j in any small 
time interval of length h is Xijh + o{h) as h ^ while the 
probability that two or more infections travel in the same 
interval is o{h) as h ^ 0); any infection arriving at node 
j changes the state of any healthy plant there to infected, 
and otherwise has no effect. All lifetimes and infection 
processes are considered independent. The initial state 
of the system is typically defined by one or more nodes 
being occupied by infected plants, the remaining nodes 
being occupied by healthy plants. The epidemic may die 
out at some finite time at which the set of infected nodes 
first becomes empty, or, on an infinite graph only, it is 
possible that it may continue forever. 

Thus, for any infected node i, the event Eij that any 
neighbouring node j receives at least one infection from 
node i has probability pij = 1 — E[exp(— A^Ti)] (where 
E denotes expectation). Note that, for any given node 
i, even though the infection processes are independent, 
the events Eij are themselves independent if and only if 
the random lifetime is a constant. We now suppose 
that this is the case and that furthermore, for all ordered 
pairs {i,j) of neighbours, we have pij = p for some prob- 
ability p. Suppose further that it is possible to observe 
neither the time evolution of the epidemic nor the edges 
of the graph by which infections travel, but only the ini- 
tially infected set of nodes and the set of nodes which are 
at some time infected and thus ultimately in the empty 
state R. It is then not difficult to see, and is indeed well 
known [TgI . [TtI , that the epidemic may be probabilisti- 
cally realised as an unoriented bond percolation process 
on the graph in which each edge is independently open 
with probability p (or closed with probability 1 — p), and 
in which the set of nodes which are ever infected con- 
sists of those nodes reachable along open paths (chains 
of open edges) from those initially infected. (Note that 
the ability to use an unoriented bond percolation process 
requires both the assumptions that the above events Eij 
are independent and that pij = pji for all i,j; in the ab- 
sence of either of these assumptions one would in general 
need to consider an oriented process with the appropriate 
dependence structure psj.) 

In the present paper we consider the epidemic to take 
place on some subset P of Z'*, where we allow P = 1^ 
as a possibility. Two sites (nodes) are considered neigh- 
bours if and only if they are distance 1 apart and thus 
Z"* induces a graph, called the d- dimensional cubic lat- 
tice h'^ = (Z'^,E'^), with the set E'' of edges connecting 
neighbours. Thus in the case P = I?, for example, each 
node has 4 neighbours. Notice that the set P also in- 
duces a graph with respect to the lattice L'^ — we shall 
denote this graph 11. This may be considered as a model 
for nearest-neighbour interaction. 
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We assume furthermore that initially there is a single 
infected site, and that all other sites in 11 are occupied 
by healthy individuals. 



B. Incomplete observations 

The probability p introduced above is considered to be 
unknown, possibly depending, through the experimental 
design, on other parameters (for example, on the spacing 
of the grid). We consider inference under each of the 
following two scenarios: 

(i) the observations consist of the set of sites which are 
ever infected, so that the routes by which infections 
travel are not observed; note that, in terms of the 
bond percolation process, this corresponds to one's 
knowledge of the connected component containing 
the initially infected site — the location of this site 
within the component not being relevant to infer- 
ence for p (see below) ; 

(ii) all that is observed is the size of the set of sites 
which are ever infected. 

We refer further to the former of these two scenarios as 
SI and to the latter scenario as 52. 

Since for any given set of initially infected sites, the 
distribution of the set of ever-infected sites in an SIR 
epidemic with constant infectious times is the same as for 
the corresponding unoriented bond percolation process 
(with the same initial set of the nodes involved), a final 
snapshot of such epidemic can be seen as an open cluster 
of the corresponding percolation process. 

Figure [T] depicts an open cluster obtained by simula- 
tion of percolation process on the integer lattice I? for 
p = 0.478. This connected component containing the 
origin can be seen as a final (and finite) outbreak of an 
SIR epidemic process of the kind discussed above. The 
origin (or, indeed, any other vertex of the open cluster) 
may be considered to be the site where the initially inoc- 
ulated individual has been placed. Clearly, the realised 
bond structure is not the only possible way which re- 
sults in the site configuration from Figure [TJ However, 
the distribution of this site configuration as an extinct 
SIR epidemic coincides with that of the corresponding 
unoriented bond percolation process. 



III. PARAMETER ESTIMATION 

A. Scenario 51 

Let n be a (proper or improper) subgraph of L"^ which 
contains the origin and let C be an open cluster of a 
percolation process on 11 containing the origin. The set 
of nodes C represents a snapshot of an extinct outbreak 
of our spatial SIR epidemic evolved on 11 C . 




FIG. 1: An open cluster (black solid dots) containing the 
origin (a black dot in a circle) as a result of percolation 
simulation on L'^. Here the bond percolation probability p 
was taken to be 0.478; the solid segments represent open 
bonds. The open cluster can be seen as a finite outbreak of 
an epidemic with constant infectious periods and neighbour- 
to-neighbour infection intensity spread rate Xij « 0.65 (since 
0.478 = 1 - e"*^ *^^). The dotted lines depict directions along 
which infection did not spread (from black to grey dots); thus, 
grey dots depict individuals which remain healthy and the 
dotted lines represent those bonds that must be closed given 
the knowledge of the cluster set. 



Let us introduce some additional notions. Let G = 
(V, E) be a locally finite graph and let G' ^ {V , E') be a 
subgraph of G (that iaV and E' C E). By the sat- 
uration of the graph G' with respect to G we understand 
the graph G = {V, E) such that 

V imAE= {{x, y)\x,y<^V'k (x, y) e E}. 

We denote the saturation of G" with respect to G by 
Saturc G' or, in cases when it is clear from the context 
with respect to what graph the saturation takes place, 
by SaturG'. A graph G' whose saturation coincides with 
itself is called a fully saturated graph. For example, the 
fully saturated graph (with respect to L^) is obtained 
from the graph depicted in Figure [T] by connecting all 
pairs of neighbouring black sites (according to the 4- 
neighbourhood relationship). Note that the operation 
of saturation may also be applied solely to a subset of 
vertices of the original graph, since it does not make use 
of the edges of the subgraph-operand. 

To distinguish between the boundary points of a graph 
and their neighbours, which are not in the graph, we 
introduce the notions of the surface and the frontier of 
the graph (again, with respect to another graph). Let us 
denote by dG the surface of G in 11, G C H, that is to 
say the set 

dG {x £ G : 3y £ n\G s.t. x and y are neighbours in H}, 

and by Fg the frontier of G in H, i.e. the set d{Il\G). 

In order to identify the likelihood function we intro- 
duce the set GiC) of all connected subgraphs of 11 with 
C as a vertex set. Note that the set Q{C) is necessarily 
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nonempty. For each G G G{C) the number of edges be- 
tween the vertices of the graph G and the elements of its 
frontier Tq is constant — we denote it by wc- Finahy, we 
denote the total number of edges in G by e(G). 

The probability that C represents the set of ever- 
infected sites and that the edges of G correspond to those 
routes along which the infection travelled is 

Pp(G') -p)'=('5aturC)-e(G)+«,c^ 

and the likelihood function associated with the observed 
set C of ever-infected sites is given by 

GeQ{C) 

Hence, under assumption of a uniform prior for p, its 
posterior distribution 7r(p | C) is a mixture of beta distri- 
butions: 

■n{p I C) oc ^ r{k) Beta (fc 1, e(SaturC) - fc + wc + 1) , 
k 

where 

r{k) ■.= #{Geg{C)\e{G)^k}. 

It is not feasible to calculate 7r(p | C) in the above form, 
since it is hard to enumerate all corresponding graphs 
and, thus, to calculate efhciently the coefficients r{k). 
We describe therefore an McMC algorithm that allows 
one to sample from the distribution 7r(p | C) under the 
uniform prior on p, that is, effectively, to estimate the 
likelihood function of p. 

Our Markov chain explores the joint space (0, 1) x Q{C) 
of possible values for p and graphs from Q{C). The sta- 
tionary distribution of the chain is the joint posterior dis- 
tribution of p and G S G{C). The description of the chain 
is given in Algorithm [T] f Appendix lA This Markov 
chain explores the set of all connected graphs Q{C) by 
simply deleting or adding an edge from the current graph 
preserving the connectivity of the given site configura- 
tion C. 

We apply Algorithm [T] to the site configuration C from 
Figure [T] (black dots only) . Figure [2] depicts the likeli- 
hood function of the model parameter for the complete 
observations (i.e. nodes and edges of the cluster) and a 
histogram of the sample obtained from the posterior dis- 
tribution 7r(p I C) when the prior distribution is uniform 
on the interval (0, 1). 

Mixing properties of Markov chains described by Al- 
gorithm [T] are briefiy discussed in Appendix |B] 



B. Scenario S2 

Under this scenario only the size n of the outbreak of 
our SIR epidemic evolving on 11 = L'^ is given. 
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FIG. 2: The solid line corresponds to the likelihood function 
evaluated for the complete information (both the sites and 
edges are known) on the cluster C from Figure [1] The his- 
togram is based on a large sample drawn from the McMC 
applied to the site configuration C (nodes only). 

Let Gn be the set of all possible connected graphs on n 
vertices including the origin. These graphs represent the 
outbreaks of the size n and we distinguish all isomorphic 
graphs which have different locations or orientations. 

We denote the number of edges of 11 between the ver- 
tices of the graph G G Qn and the vertices of its frontier 
Fg by w{G). 

Given the epidemic size n, the inference on p involves 
evaluation of the likelihood function Cn{p) ■= Fp{\C\ = 
n) which can be represented as follows: 

As previously, under assumption of a uniform prior for 
p its posterior distribution n{p\ \C\ ~ n) is a mixture of 
beta distributions: 

7r(p I |C| = n) oc ^ q{s, fc, I) Beta {k + I, s - k + I + I) , 

s^k.l 

(1) 

where 

g(s, k, I) := #{G e Gn I e(Satur G) = s, e(G) = k, w{G) = /}. 

This, again, represents a hard enumeration problem. 
However, inference on p can be made using the McMC 
technique. Algorithm [2] (Appendix lA ip contains a de- 
scription of a Markov chain which serves the purpose of 
sampling from the posterior distribution 7r(p| |C| — n), 
given the prior distribution of p is uniform on [0,1]. The 
chain explores the joint space (0,1) x Gn of possible val- 
ues for the percolation parameter p and all possible con- 
nected graphs on n nodes, by deleting a vertex from a 
graph and adding a vertex from its frontier. 
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FIG. 3: Likelihood functions Cn(p) {n = 25,50,70) obtained 
using the McMC from Algorithm [2] and McMC sample his- 
togram of £n{p) for n — 70. 



Figure [3] depicts plots of the likelihood function 
when n = 25, 50, 70 and a histogram of a sample drawn 
from tt{p \ \C\ — 70). The graphs of the likelihoods were 
obtained by smoothing the histograms of correspond- 
ing samples generated by the McMC presented in Al- 
gorithmic] Note that, as might reasonably be expected, 
the values of p maximising the likelihood increase with 



Mixing properties of Markov chains described by Algo- 
rithm [5] are briefly discussed in Appendix [Bl Two video 
examples of dynamic graph updates illustrating Algo- 
rithms [1] and [2] can be found at the WEB address [18]. 



C. Convergence analytical results for <S2 



1. Asymptotic behaviour of the maximum likelihood 
estimates of p 

Denote by C{x) the open cluster (connected compo- 
nent) which contains the vertex x. Let us write xip) — 
Ep|C| for the mean number of vertices in the open clus- 
ter at the origin. Using the translation invariance of the 
process on Z"*, we have x(p) — Ep|C(a;)| for all vertices x. 
Percolation theory tells us that if p < pc, then x(p) < oo. 
When p > pc, then xip) = oo and the function x is 
not of a much interest in this case. Instead, one studies 
the function x^ [p) = ^p[\C\ ■ \C\ < oo]. The function 
x{p) iX'^ {p)) monotonically increases to infinity as p t Pc 
{p -l Pc)- It is known [Toj that there is no infinite open 
cluster when p = pc for percolation on the square lattice 
and on the lattice Z'', where d > 19; this is also be- 
lieved to be true for all other values of d. How likely it is 
to observe an open cluster of size n when n is very large? 
What value of p should one suggest if one happened to 
observe a large epidemic of size n? 

Intuitively, it is unlikely that very large but finite epi- 
demics may be obtained for values of p significantly dif- 
ferent from Pc (if p Pc, then it is unlikely that the 
epidemic will attain a very large size; if p ^ Pc, then it is 
equally unlikely that, having attained a sufficiently large 
size n, the epidemic will have burned out). Intuition sug- 
gests therefore that the likelihood function for p, given 
that the size n of the connected component containing 
the origin is growing, should be increasingly concentrated 
around Pc- 

Let Cn{p) = IPp(|C| = n) be the probability that an 
open cluster C from a percolation process with the edge 
density p is of size n. (Equivalently, assume that we 
observe an SIR spread of infection via nearest-neighbour 
interactions on and that the final size of its outbreak 
is n; when n is fixed the function Cn{p) can be regarded 
as the likelihood function for the percolation probability 
parameter p.) 

Let Pn be the maximum likelihood estimate for p 
derived from £„(p). For any p define also Ln{p) — 

Cn{p)/ Cn{Pc)- 

Theorem III.l The sequence of maximum likelihood es- 
timates Pn for p converges to the critical probability pc- 



Percolation exhibits a phenomenon of criticality, this 
being central the percolation theory: as p increases, the 
sizes of clusters (connected components) also increase, 
and there is a critical value of Pc at which there appears 
a cluster which dominates the rest of the pattern, li p < 
Pc, then with probability one all open clusters are finite, 
but when p > Pc there is a single infinite open cluster 
almost surely Note, that the critical percolation 

parameter pc depends on the dimension d of the lattice. 
Bond percolation on the square lattice seems to be most 
studied to date of all percolation processes. The critical 
probability Pc in the case of a square lattice is ^ . What 
follows, however, holds for any lattice Z'^, d > 3. 



The proof of this theorem is based on the following 
lemma, together with a generally accepted regularity con- 
dition. 

Lemma III. 2 For any p G (0,1) different from Pc the 
following holds: 

liminf Ln{p) — 0. 

n~¥oo 

Moreover, this convergence is uniform for any closed in- 
terval which does not contain pc- 

The proof of the lemma and the theorem can be found 
in Appendix I A 21 
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As mentioned above, it is reasonable to expect the se- 
quence of the maximum hkehhood estimates p„ to be 
monotonic (see Figure |3|). This observation gives rise to 
the fohowing conjecture. 

Conjecture III. 3 The sequence o/{p„} converges to Pc 
monotonically from the left. 

Furthermore, we formulate the corresponding compan- 
ion conjecture for the sequence of posterior distribution, 
using the notion of a delta sequence . 

Conjecture III. 4 Provided p £ suppn^-) the functional 
sequence {7r(p | n)}^^]^ is a delta sequence which gener- 
ates the delta function 5{p — Pc)- 

Thus, we believe that the limiting posterior distribution 
of the percolation parameter is a one-point mass distri- 
bution at p = Pc, or the Dirac delta function S{p — pc). 

2. Combinatorial characterisation of large percolation 
clusters on h"^ 



that is the number of open edges in G is approximately 
equal to the total number of closed edges and edges be- 
tween G and its frontier. 



3. Large percolation clusters as rare events 

When n is large, the appearance of finite open clusters 
of size n is highly unlikely: the distribution of the clus- 
ter size (hypothetically) decays as n"^"^/*, S > 0, when 
p = Pc, and the decay is exponential (sub-exponential) 
when p < Pc {p > Pc). Large finite percolation clusters 
can therefore be viewed as rare events. Since the state 
space of the McMC proposed for inference on the per- 
colation parameter p under scenario S2 and described 
in Algorithm [2] involves the set of all open clusters on 
n nodes, this algorithm can be readily used in order to 
obtain realisations of these rare events. 



IV. BAYESIAN OPTIMAL DESIGNS 



The theoretical results obtained and conjectured pre- 
viously for inference under scenario 52 can be used to 
derive their combinatorial analogues regarding the rela- 
tive number of realisations of the process with the clus- 
ter size n. Under scenario 52 the posterior distribution 
7t{p I \C\ — n) can be seen as a mixture of beta distribu- 
tions, as in ([Ij. As in Section fill Bl let q{s,k,l) be the 
number of all possible graphs G which may correspond 
to an open cluster C (the vertex set of G and C coincide) , 
and such that k is the number of edges in G, s is the 
number of edges in the saturation of G, and I is the num- 
ber of edges of H between the surface and the frontier 
of G. Conjecture IIII.4I would suggest that the number 
q(s, k, I) of graphs G corresponding to open clusters C 
which could emerge as a result of the percolation process 
with parameter pc and for which it holds that 



k + l 
s + l + 2 



Pc, 



(2) 



is far greater than the number of all other graphs. This is 
so, since the sequence of beta distributions Beta(a„,/3„) 
is a delta sequence generating the delta function at pc if 
and only if a„/ (a„ + Pn) Pc- 

Thus, in percolation processes on L'' the number of 
finite graphs corresponding to open clusters of size n 
(where n is large) that satisfy the condition 



e(G) + 1 



e(SaturG) + w(G) 



Pc 



(3) 



largely exceeds the number of all other connected com- 
ponents on n nodes. In other words, a typical graph 
corresponding to an open cluster of a large size in per- 
colation process on L'' is a connected subgraph of L'' 
characterised by ([3]). In particular, when d = 2: 



A. Utility based optimal designs 

1. General description 

The experimental design problem can be conveniently 
approached within the Bayesian framework 13] . Suppose 
we study a stochastic process for which we formulate a 
model M, characterised by a model parameter 6. The 
model M specifies a sampling density fd{y \ 0) for the ex- 
perimental outcome y obtained from observation of the 
process under experimental set-up d given the value of 
the model parameter 0. Our knowledge about is de- 
scribed by a prior distribution tt{6). Whenever the choice 
of set-up d is under our control there arises the question 
of selecting the optimal d under which one should observe 
the stochastic process. Such prescribed conditions are re- 
ferred to as a design, and the optimal design is found un- 
der optimality criteria specifically tailored to the context 
and the purpose of the experiment. 

By employing a utility function u{d, y, 9) one can spec- 
ify the purpose of the experiment and measure the value 
of its outcome y accordingly. The methodology of pos- 
ing and solving utility-based optimal design problems 
within the Bayesian paradigm became somewhat stan- 
dard [3 The design has to be chosen before per- 
forming the experiment and one may choose to maximise 
the expectation of the utility function u{d, y, 0) with re- 
spect to 9 and y [l^ : 



where 



e(G) - w{G) w e(Satur G) - e(G); 



(4) 



U{d) 



argmax{7(c?), (5) 



u{d,v,9)fd{y\9)n{9)d9dy. (6) 



e Y 
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Here V is the set of possible designs. The set of possi- 
ble outcomes y of the experiment is denoted by Y . The 
experiment is defined by a model fd{y \ 0), i.e. by the sam- 
pling density of y conditional on 6 for a given design d. 

When the purpose of the experiment is to infer the 
model parameter, a sensible choice to measure the utility 
of the experiment under design d might be one that rep- 
resents the information gained on 9 from the experiment. 
One of the standard choices is the KuUback-Leibler (KL) 
divergence £>kl between the prior distribution of 0, 7r(0), 
and its posterior distribution Tr{0 \ y, d) obtained follow- 
ing the experiment yielding data y: 



DKL{7Ti0\y,d) II 7r{e)} 



io,^-^AO\y,d)d0, 



(7) 

where 71(91 y,d) (X fdiy\9)n{0). 

Introduced by KuUback and Leibler [22| this in- 
formation measure was justified by Lindley [231 and 
Bernardo [24] and is related to the Shannon information 
measure. 

The KL divergence is a random variable, being a func- 
tion of the random outcome y. In order to find the the 
design dmax that is optimal with respect to this choice of 
utility measure one should maximise the expected value 
of the KL divergence over the set of observables y as a 
function of the design variate d: 



UKL{d)=^Y[DKM0\y,d) II 7t{0)}], 



(8) 



this being invariant under a reparametrisation of the 
model in terms of 0. 

When the experimental motivation consists in increas- 
ing one's knowledge on the model parameter the ex- 
pected KL divergence UKhid) coincides with ^ taking 
u{d,y,0) = log ^^^^^^jy^ . If, however, the purpose of the 
experiment is to inform someone who holds a different 
prior to our own and we wish to advise them on which 
design to use, using our own 'superior' knowledge of the 
system under study, then the expected KL divergence 
should be calculated in the form prescribed by (|5]) when 
the integration over the space of observables Y is to be 
carried out using one's superior knowledge. We refer to 
these contrasting experimental scenarios as progressive 
design and instructive design [2^, respectively, and dis- 
cuss them next in more detail. 



2. Progressive design 

Under this scenario there is an experimenter who 
holds prior beliefs regarding the parameter, 0, in the form 
of a prior distribution n(0) (perhaps obtained from earlier 
observation of the process) and whose goal is to design an 
optimal experiment in order to increase this knowledge. 
Using KL divergence as the measure of information gain, 
this experimenter should maximise the expected KL di- 
vergence ([S]): 



UKL{d) = EyEc- 



lo, 



7Ti9\y,d) 



log 



e Y 



AO) 

<0\y,d) 



fd{y,9) d0dy, 



so that UKhid) coincides with © where u{d,y,0) = 
log ^''^(g^^ ■ Notice, however, that the expected informa- 
tion gain UKhid) can be written as follows: 



UKhid) 



Y e 



log7Ti0\ y,d)fdiy\0)7ri9) d0dy 



log7r(6')7r(6')d6l. 



where the second term is the Shannon entropy of 7r(0), 
which does not depend on d. Thus, in order to obtain 
the solution cJmax in this case, it suffices to maximise the 
first term only. 



3. Instructive design 

In contrast to the progressive design scenario, in the in- 
structive case there is an experimenter S", holding a prior 
7r(0), and a better (or differently) informed trainer (or 
instructor) £^ whose knowledge about the model param- 
eter is summarised in a distribution n* (0). (Of particular 
interest is a special case where the instructor knows the 
true value of 0.) The aim here is to maximise the change 
in experimenter's belief from 7r(0) to 7r(0 | y) by designing 
an experiment using the superior knowledge 7r*(0). 

This latter optimisation problem can be formulated as 
follows: 



Cax = argrnaxC/j^Llc')- 



(9) 



C/KL(rf) - / DKhM0\y,d) II 7Ti0)}riy)dy, (10) 



where, as before, 77(6* [ y) is derived as in d?]), and f*iy) 
is the instructor's prior predictive distribution for y: 

riv)^ J fdiy\9)n*i9)d9. 
e 

In particular, if the instructor £^ knows the exact value 
of 0, 9*, and hence 7r*(6') is the Dirac function 6i0 ~ 0*), 
then r(y) = /d(y|r), so that 

UKhid) - I DKh{ni0 I y, d) \\ ^(^)}/rf(y | 9*)dy. (11) 
y 

Note that the 'instructive' expected utility C/^l gen- 
erally not symmetric with respect to the exchange of 
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tt{9) and tt*{9): if the experimenter (o and the instruc- 
tor £^ change their roles, then the corresponding optimal 
instructive designs are different, unless their knowledge 
about 9 is identical, that is Tr{9) ~ 7r*(f?) almost every- 
where. 



B. Simulation-based evaluation of the expected 
utility 

Generally, the solution to the optimal design problems 
(l5]|6|) and (|9l lT0l) cannot be obtained analytically. This 
is mainly due to the following three problems. First, 
the design space V can be complicated, with many de- 
sign variables, some of them having a continuous range 
of values. Second, even if the design space has a sim- 
ple structure, the utility function may not be simple to 
evaluate, so the expected utility U{d) cannot be obtained 
explicitly. Finally, for incomplete observations of highly 
non-linear stochastic processes such as epidemic models, 
the likelihood is not usually available in a closed form, 
and this results in computationally intensive evaluations 
or sampling methods. 

A review of analytical and approximate numerical solu- 
tions to Bayesian optimal design problems for traditional 
experimental design involving linear and non-linear mod- 
els can be found in [l^, [l^ . 

Miiller [3| reviews simulation-based methods for op- 
timal design problems where the expected utility U {d) 
is evaluated by Monte-Carlo simulation. In its simplest 
form an estimate U oi U for any given design d in the 
progressive case is as follows: 

1 

where {{9i,yi),i — 1,...,M} is a Monte-Carlo sample 
generated values: 

9^,^7r{9), y,^p4y\9). (13) 

The expected utility U may, in particular, be based on 
the KL divergence. 

Analogously, the expected utility C/j^L under the in- 
structive scenario, when the instructor knows the true 
value of the model parameter, 9*, can be evaluated using 
the following scheme: 



1=1 

where y* - f{y\9*,d), i = 1,...,M, and KL{y^,d) 
is an estimate of the KL divergence DKh{T^{9 \y* , d) \\ 
7r(^)} that can be obtained via numerical integration 

of log ^^^^1^^ with respect to the posterior 'K{9\y,d). 
The former function, in turn, may need to be evaluated 



through simulation methods — perhaps using an McMC 
scheme. 

Evaluation of a continuous expected utility surface 
U{d) by or (|14p can be achieved by computing its 
values on a discretised grid of points and further smooth- 
ing of the obtained set of values in order to approximate 
the expected utility landscape. This, however, may be 
problematic when d e 2? is a multidimensional design pa- 
rameter. An alternative method is provided by the aug- 
mented probability simulation approach which is studied 
in m, [2^ , and reviewed in [3] . 

The augmented probability simulation approach as- 
sumes that M(d, 6*, y) is a non-negative bounded function. 
(This condition, although not always automatically satis- 
fied, can be easily achieved by correspondingly modifying 
the utility function.) An artificial density, proportional 
to u{d,9^)fd{y\9)TT{9) can then be defined for {d,9,y) 
jointly Q: 

h{d,9,y)^uid,9,y)Uy\9)7r{9), (15) 

so that its marginal in d is proportional to the expected 
utility U{d). Sampling from h{-, ■, ■) can be used to obtain 
a sample from its marginal using the Metropolis-Hastings 
McMC scheme described in [l4|. 

When using the augmented probability modelling ap- 
proach one approximates the optimal design by an empir- 
ical mode of the marginal density of d under the artificial 
distribution h, that is by the mode of histogram formed 
from the first components of the set of samples obtained. 
Notice that the evaluation of M(d, 6*, y) at each sampling 
step may in general itself involve an McMC sampling 
from the posterior 7r(0 \ y^d)^ so that the approach is po- 
tentially very expensive computationally. 

C. Inner-outer plots 

We now apply the utility-based Bayesian approach in- 
troduced above in order to design optimal experiments 
for our SIR epidemic model with incomplete observa- 
tions of the kind described by scenario SI. We do so 
by identifying a (finite) class of designs (node configura- 
tions) which we call inner-outer plots. These are plots 
of bounded size which are obtained by removing some 
nodes from the underlying grid. We refer to this process 
as sparsification of the grid. 

A typical example of an inner-outer plot in 1? is given 
in Figure m The inner-outer plot depicted consists of two 
parts — inner and outer regions (hence the name of the 
design). Some sites are removed from the outer region 
of the plot in a regular manner, while all sites of the in- 
ner region are preserved. Our motivation in considering 
sparsified grids is to show that intuitively obvious de- 
signs, such as completely dense grids, can be improved. 
The justification of the grid sparsification is that it im- 
pedes the ability of the epidemic to spread, which makes 
possible better inference when p is large, while the more 
dense inner part of the design works better when p is 
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FIG. 4: An example of an inner-outer (m, r)-plot (nodes only) 
in Z^. In this inner-outer plot m = 9 and r = 3. The plot is 
bounded by an A*' x A'' square, where A^ = 21 can be calculated 
using (lisp . Depicted is the saturation of an (9, 3)-plot with 
respect to L^. 



small. A formal description of inner-outer plots follows 
next. 

Let us assume that m is an odd positive integer and 
r G N U {0}. An inner-outer {m,r)-plot IIq'*'' (m,r) in 

TL'^ with centre at the origin is a c?-dimensional box i?^^ 
with side-length N = m + Ar and some vertices removed 
as follows: 

Bn \C}^ {'>TT-^r), r>0, 
where C^' (to, r) is the set 

{x e B^^^ : II X Hoc- TO + 2j + 1, j = 0, . . . , r - 1 & 
& II X ||i=0mod2}. 

d 

Here || x ||i= and || x ||oo= max(|a;i|, . . . , \xd\) for 

i=l 

any x =(xl^ ■ ■ ■ , Xd) G and the box -B^'' is defined as 
follows [201 : 

: = [-{N - l)/2, [N - 1)/2Y 
^{xeZ'' ■.\\x\\oo<{N -l)/2}. 

We call any plot that can be obtained by translating the 
plot Ilg''' (m, r) in an inner-outer (to, r)-plot, or simply 
an inner-outer plot, and denote it by Il^'^^m, r). 

The total number of nodes contained in an (to, r)-plot 
can be calculated by subtracting the total number of the 



nodes removed from the outer plot (see ((16])) as follows: 

r = 7V(TO,rf-4^(^+2z-l) 
1=1 ^ ^ 
= {m + 3r)^ +r{3r + 2), (17) 

where, as before, 

A^(TO,r) = TO-|-4r. (18) 

The inner-outer (to, r)-plot depicted in Figure|3]is from 
Z"". The size of the inner plot is to x to, where to = 9, 
and there are r = 3 'circles' (with respect to the metric 
II • ||oo) in the outer plot from which every second node is 
removed. The size of the bounding box is N x where 
N = m -\- Ar = 21. The total number of nodes T that 
this node configuration contains can be calculated using 
(IT71) and is equal to 357. 

D. Inference for SIR epidemics on inner-outer plots 

We consider an SIR epidemic with constant infectious 
times which arises at the central site of an inner-outer 
plot n = n''') (to, r) from Z'^ and evolves on that plot 
according to the nearest-neighbour interaction rule, that 
is using edges from L'* with endpoints from H, subject 
to being restricted by the plot boundary. We consider 
this model in the context of scenario 51, that is when 
the only information available about the outcome of the 
epidemic is its site configuration C. As before, this model 
is equivalent to that of percolation on n^'^)(TO, r). 

We note that Algorithm [T] and Algorithm [5] can be 
readily used for making inference on the percolation 
probability parameter p given the sites of an open cluster 
C, or merely its size in the case of Algorithm[2l generated 
from a percolation process on any locally finite graph [3l| . 

Figure [5] shows the plot of a simulated open percola- 
tion cluster C obtained on the inner-outer plot H^^^ (13, 2) 
with p = 0.9 (left plot) and an estimate of the likelihood 
formed from a histogram estimate of the posterior den- 
sity of p obtained using the Markov chain described in 
Algorithm [T] with C as the input data, with a uniform 
prior density for p. Figure [5] shows similar plots for a 
realisation of the percolation process on n(^)(23,4) with 
p = 0.86. 

E. Optimal design problem and design space 

We consider now an optimal design problem for the 
percolation process on node sets of limited size. We adopt 
the utility-based Bayesian approach formulated in Sec- 
tion |IV3 The choice of the design space V can be made 
in a number of ways, possibly refiecting such restrictions 
as limitations on the number of experimental units or 
the dimensions of the experimental plot. In the context 
of inner-outer plots the former restriction would mean 



10 




FIG. 5: Left: site configuration obtained on the inner-outer (13, 2)-plot with p = 0.9. Right; McMC histogram of an McMC 
sample for p assuming a uniform prior for this parameter. 
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FIG. 6: Left; site configuration obtained on the inner-outer (23, 4)-plot with p = 0.86. Right; histogram of an McMC sample 
for p assuming a uniform prior for this parameter. 



that T from (fT7|) is bounded, whereas the latter condi- 
tion is equivalent to bounding the quantity N{m,r). A 
combination of these conditions or some other informa- 
tion can be also taken into account when identifying the 
design space. One advantage of grid-based plot designs 
is that they can greatly reduce the dimensionality of the 
design space compared with designs that allow units to 
be placed at arbitrary locations. 

We assume now that N is odd and fixed and that the 
design space has the form 

V= {d\d = n(2) (to, r) & m + 4r = TV} 

For example, if iV = 19 then, as can be derived easily 
from (fT5)) . the design space consists of the following de- 
signs: 

p = {n(2)(i9,o), n(2)(i5,i), n(2)(ii,2), 
n(2)(7,3), n(2)(3,4)}. 

The set of observables Y is the set of all connected com- 
ponents on d g D containing the central node. 

V. PRACTICAL IMPLEMENTATION 

Given a finite design space T) it is fairly straightfor- 
ward to solve the optimal design problem for the percola- 



tion model on inner-outer plots using the tools developed 
above. We now make some comments about solving the 
problem under each of the two design scenarios. 



A. Progressive design: expected utility evaluation 
through augmented modelling 

Since V is finite we choose to identify the design that 
maximises the expected utility function in the progres- 
sive case using augmented modelling as described in Sec- 
tion |lVBl Recall that this is based on an artificial dis- 
tribution h{d,p,y) oc u{d,p,y)fd{y\p)Tr(j)) from which 
samples are generated using a Metropolis-Hastings sam- 
pler. The optimal design d* is identified then as a value 
of d at which the marginal of h is maximised. 



B. Instructive design: Monte Carlo evaluation of 
the expected utility 

We treat the 'instructive' case differently from that of 
the 'progressive' one because of the form of the expected 
utility under this scenario. Recall that in the 'instructive' 
case whenever the instructor knows the true value p* of 
the model parameter p one can write the expected util- 
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ity based on the Kullback-Leibler divergence as follows 
(compare with (ITT]) '): 

1 

C^klW = E -MylP*) I log ^^^^<P I y, d) dp, 

where fd(y\p*) is the likelihood function evaluated at 
the true value of the model parameter p* given the open 
cluster y at the centre of the inner-outer plot d. 

It is clear that to evaluate the expected utility C/j^L (^) 
via standard Monte Carlo simulation, then a Markov 
chain has to be run to evaluate the posterior TT{p\y,d) 
each time we sample for a new observation (an open 
cluster) y and also the potentially time-consuming in- 
tegration has to be done with respect to the model pa- 
rameter p. This integration can be implemented in the 
following way. Since we can sample from the posterior 
7T{p\y,d) via McMC (Algorithm [TJ for any given open 
cluster we do so and then fit the beta distribution (or 
some other distribution) to the McMC sample obtained 
in order to perform integration numerically in a more 
efficient , albeit, approximate manner. 

Thus, the expected utility evaluation scheme for inner- 
outer plots in the 'instructive' case and scenario 51 can 
be described as follows. 

For each inner-outer plot d G I? we do the following: 

(i) generate a random sample of M independent 
connected clusters {yi}f£i on U{m,r): yi ^ 

f{y\p*,d); 

(ii) perform AI McMC's in order to obtain AI in- 
dependent samples for the posterior distribution 

(iii) fit a Beta distribution to each of the sam- 
ples obtained; refer to the fitted distributions as 
■^{p\yi,d); 

(iv) evaluate numerically the integrals 



I. ■■= / log 



.T^{p\yi,d) 



n{p) 



7r{p\y,,d)dp] 



(v) estimate the expected utility: Um 



M 



It is natural to ask how well the true posteriors can 
be fitted by a Beta distribution (step (iii) of the above 
scheme). Experience suggests that such a fit never affects 
the outcome of the analysis on the qualitative level unless 
the prior distribution has more than one local mode or 
its support is smaller than the entire interval (0, 1). Re- 
call that the purpose of this step is to make evaluation 
of the integrals Ii easier and faster, and hence the family 
of Beta distributions is only one possible choice. For in- 
stance, if the prior distribution 7r(p) is uniform on (0, 1), 
then the integrals Ii represent the entropies of the fitted 
distributions. In the case when the fitting is done by 
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FIG. 7: Inner-outer design plots A, B, and C form the design 
space T> = {A,B,C}. 



beta distributions, each of these integrals can be quickly 
calculated using the following analytical formula for the 
entropy of the beta distribution Beta(Q;, /3): 

Ent{Beta(a, (3)} = log B{a, 13) + {a + l3 ~ 2)V'(a -t- /3) 
-(a-l)V'(«)-(/3-l)V(/3), 

where tp is the digamma function, i^{z) = F'(a;)/r(2:). 
Other choices of prior distribution may suggest alterna- 
tive and more appropriate fitting distributions to facili- 
tate the calculation of the integrals Ii, i — 1, . . . , M. 



C. Example 

In our example we consider all inner-outer plots in 
whose sizes do not exceed = 11. There are only three 
such plots: n(2)(3,2), n(2)(7,l), and n(2)(ll,0). For 
ease of reference we mark them A, B, and C respectively 
(as depicted in Figure [7]). Thus, the design space V = 
{A, B, C} consists of three designs, among which A is the 
most sparsified plot whereas no nodes are removed from 
C at all. 

Figure |S] represents graphically the results of the com- 
parison of designs from 2? under both 'progressive' and 
'instructive' scenarios when the prior distribution n(p) 
is uniform on the interval (0,1). The left panel of the 
figure corresponds to the former scenario and depicts a 
histogram of a sample corresponding to the marginal 
of the artificial augmenting distribution h{d,p,y) cx 
u{d,p,y)f{y\p,d)Tr{p) in d € T>. The right panel corre- 
sponds to the latter scenario and shows the Monte Carlo 
estimated values of the expected utilities and 95% cred- 
ibility intervals for each of the three considered designs 
(Af = 1500, see ([H]) in Section ITVB)) assuming that the 
instructor knows p precisely so that 7r*(p) = S{p — 0.9). 

The plots from Figure [5] indicate that the solutions 
to the optimal design problem under the two scenarios 
are different from each other. The 'moderately sparsi- 
fied' plot B maximises the expected utility in the pro- 
gressive case, that is in the case when the experiment is 
designed by a single experimenter. If, however, it is the 
instructor who knows the true value of the model param- 
eter {p = 0.9) and wishes to choose the best inner-outer 
plot from the set T> for a more uninformed experimenter 
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FIG. 8: Left: sample histogram for the marginal of h{d,p,y) in d, d £ {A, B,C}, under progressive design and 7r(p) ~ U(0, 1). 
Right: evaluated expected utility under instructive design with 7r*(p) = S{p — 0.9) and 95% credibility intervals (Af = 1500) 
for the plots A, B, and C, under instructive design. 



to use (instructive scenario), then the optimal plot is 
the 'mostly sparsified' inner-outer plot A. Notably, the 
'mostly dense' plot C would be the worst choice in the 
instructive case, whereas it outperforms the 'mostly spar- 
sified' plot A but is worse than the 'moderately sparsified' 
plot B in the progressive case. These results demonstrate 
the moderately counter-intuitive finding that the optimal 
design is not necessarily that with the maximal number 
of experimental units and that a 'less-is-more' principle 
may be in operation. On another level, this finding may 
not be entirely unpredictable in the instructive case. The 
optimal design's combination of a dense centre and spar- 
sified outer region may help to increase the probability 
that the experimental epidemic succeeds in establishing 
itself but does not lead to saturation of the lattice before 
it dies out. This in turn may represent an outcome that 
yields much information on the unknown p. 



VI. CONCLUDING COMMENTS 

In this paper we have considered a number of question 
relating to the observation of percolation processes on fi- 
nite lattices. We have demonstrated the value of McMC 
techniques for parameter estimation under different sam- 
pling regimes and shown how Bayesian experimental de- 
sign techniques can be used to identify the optimal lat- 
tices on which to carry out experiments. 

We introduced inner-outer plot designs by sparsifying 
the underlying grid. Allowing one to control excessive 
connectivity in outer regions of experimental plots, these 
designs help to separate different values of the model 
parameter p when p ^ pc, while giving the epidemic 
a chance to establish in inner regions when p ^ pc- 
Inner-outer plots can be further generalised by consider- 
ing more than two levels of grid sparsification, according 
to the principle used in our work: the farther the region is 
from the centre of the plot the more intensive sparsifica- 
tion it receives. To an extent, this process can be viewed 



as switching to similarly shaped lattices with percolation 
parameter taking its values from the set {p,p^,p'^, . . .} in 
decreasing fashion. 

Although we have worked largely within framework 
of percolation theory, the connection between percola- 
tion and spatio-temporal epidemic models with nearest- 
neighbour interactions mean that the results may have 
significant implications for and application to the study 
of epidemic systems. 
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Appendix A: Algorithms and proofs 
1. McMC algorithms for the scenarios 51 and <S2 

An McMC algorithm for making inference on the per- 
colation probability p under scenario 51 is given in Al- 
gorithm [TJ The proposed McMC is irreducible by con- 
struction: there is a positive probability for the chain 
to switch between any two connected graphs from G{C), 
since any two such graphs share a common vertex set and 
differ only by a finite number of edges. 

Algorithm [2] describes an McMC algorithm for mak- 
ing inference on p under scenario 52. This algorithm, 
similarly to Algorithm [1] is a combination of Gibbs and 
Metropolis-Hastings steps. The chain's stationary distri- 
bution f{p,G) specified a marginal density for p which 
coincides with the posterior distribution tt{p \ \C\ = n). 

We give now explicit expressions for the proposal prob- 
abilities used in the Metropolis-Hastings steps within 
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Algorithm 1 Markov Chain Monte Carlo: scenario SI Algorithm 2 Markov Chain Monte Carlo: scenario S2 



Require: an open cluster C; 
1: take an initial value po arbitrary from (0, 1 
t--OXt ■- (pt, SaturC); 
repeat 

Gibbs sampler steps: 

'+ l,e(SaturC) 
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e(G2t) + ™c + l); 



P2t+i ~ Beta(c(G2t 

X2t+1 (p2t+l,G2t)i 

Metropolis-Hastings sampler steps: 

choose an edge e uniformly at random from SaturC; 

if e e G2t then 

(7 ~ Uniform [0,1]; 

it U < min(l, ^^^^l{G2t+i\{e} is connected}) then 

X2t+2 ■■= {p2t+i,G2t+i \ {e}); 

else 

^2t+2 := (P2t+1, G2t+l); 

else 

f/ ~ Uniform [0, 1]; 

ifU< min(l, /"'+^ ) then 

— \ ' 1-P2t+1 ' 

^2t+2 := (p2t+i,G2t+i U {e}); 
else 

X2t+2 (P2t+1, G2t+l); 

t:=t + 2- 

until we judge that the chain has converged and a sample 
of sufficient size is obtained 



this algorithm. Assume that the current graph within the 
Metropolis-Hastings step is G and a graph G is proposed, 
the latter being obtained from the former by deleting a 
vertex u with all edges adjoining it and inserting a vertex 
V, including each potential edge to v independently with 
probability p, this probability having been drawn in the 
preceding Gibbs step. We assume that at least one such 
edge is inserted and denote the number of deleted and 
added edges by d{u) and d{v) respectively. Then, 



Require: the value of n. 
1: take a value po arbitrary from (0, 1) and a graph Go ar- 
bitrary from Qn; 



t~OXt := (pt,Gt); 

repeat 

move_on:=l; 

Gibbs sampler steps: 

P2t+i ~ Beta(e(G2t) + l,e(SaturG2t) - e(G2t) + w(G2t) + 1); 

X2t+l ■- {p2t+l,G2t)\ 

Metropolis-Hastings sampler steps: 

choose a vertex u uniformly at random from G2t+i and 



choose a vertex v uniformly at random from 



De- 



rive a graph G from G2t+i by deleting all edges which 
adjoin u (in G2t+i) and adding the edges that connect v 
with vertices of the graph G2t+i \ {u\ in D, each indepen- 
dently with probability P2t+i (conditioning on the event 
that at least one edge is added), 
if G is disconnected then 
^2t+2 ~ (P2t+i, G2t+i); move_on:=0 
if move_on then 
U ~ Uniform(0, 1); 

d[v) ■- #{e I e = [v, z) 3z £ G2t+i \ {u}}- 
d{u) ■- #{e I e = (w, z)3z(£G\ {v}}; 

u{u) ■- #{x I X e TG2t+i & (u, x) e n}; 

iy{v) := #{x \x €r^&i{v,x) € H}; 

K := d{u) — d{v) -\- ^{v) — v{u); U ~ Uniform(0, 1); 



19: if U < 



^ ircat + il l-(l-P2t+i)'''"' 



-(l-P2t+l)''<"' 



'P2t + l)' 



then 

20: X2t+2 ■- (P2t+1,G); 
21: else 

22: X2t+2 (P2t+1, G2t+l); 
23: t ■.= t + 2- 

24: until we judge that the chain has converged and a sample 
of sufficient size is obtained 



Hastings step, a, is as follows: 



(r r\ - 1 1 /(")(l-p)'^»-'^(") ^ 

9^^. n\TG\ 1 - (1 - p)div) connected}. 

(Al) 

and similarly, 



~ 1 1 /H(l 

q{G,G) — — T— — r -~— ]1{G is connected}- 



(A2) 



Clearly, 



Pp(G) oc/(^)(l -p)'^^('")(l 
so that the acceptance probability at the Metropolis- 



. / qiG,G) FpiG)\ 

a — mm 1, — ^ 

\^ 'g(G,G)Pp(G) J 

- mill (l I^gI (1 - p)''(")+'^(^) 1 - (1 -p)''» \ 

" ™™ V ' (1 - 1 - (1 - py(") J 

^min,fl,^i-li^(l-,)'=V 

where, as it was introduced in the description of Algo- 
rithm H 

K — d{u) — d{v) + v{v) — i'{u). 

We claim that the constructed chain is irreducible. In 
terms of graph theory this means that the chain can pro- 
ceed from each connected graph on n vertices, including 
the origin, to any other connected graph including the 
origin of the same size on the underlying lattice. We 
show the irreducibility of the proposed McMC by con- 
structing a sequence of steps in which any graph of Qn is 
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transformed to a so called line-skeleton graph on n ver- 
tices. By such a graph we mean any tree (a graph with 
no cycles) containing the origin and having n vertices so 
that only two of these vertices have degree one. Select 
one such line-skeleton and denote it by S. 

Consider a graph G from Qn and denote the length of 
the shortest open path from x G G to S hy S{x, S) (chem- 
ical distance). Each vertex x from G receives a well de- 
fined finite weight S{x, S), since the graph G is connected 
and finite. By using the description of our Markov chain 
we can delete any vertex from our current graph for which 
S(x,S) is maximal and add to this graph a vertex from 
the chosen line skeleton S without making the graph dis- 
connected or containing cycles until the maximum value 
of d{x, S) is zero — in this case all vertices are forming 
the line skeleton. Since this procedure can be reversed it 
follows that Qn in the described Markov chain is indeed a 
communicating class, and hence the chain is irreducible. 

2. Lemma and Theorem 

Proof. fLemma lIII.2[) We consider the following two 
cases [HI: 

(i) Subcritical case p < Pc- In this case the cluster size 
distribution Ln{p) decays exponentially, i.e. 

3A(p) > : Cn{p) < e-"^(P) Vn > 1. (A3) 

However, it is also known that at p = pc the 
cluster size distribution, while almost surely finite, 
nevertheless has an infinite mean. This implies 
that, for any A such that < A < Ap, we have 
lim sup e^"' Cn (pc) — oo, which when combined with 

(|A3|) . gives liminfL„(p) = 0, as required. Uni- 

n— f oo 

formity of convergence on any closed interval con- 
tained within (0,Pc) follows since A above may be 
chosen so that Ap is uniformly bounded away from 
A for p belonging to such an interval. 

(ii) Supercritical case p > pc- The argument here is 
essentially the same. In this case the decay is sub- 
exponential: 

3 77(p)>0: £„(p) < e-"'")"'""'"" Vn > 1. (A4) 

The infinite mean of the cluster size aX p = pc also 
implies that for any 77 such that < ?/ < 77^ , we have 

limsupe''"' ''' Cn{pc) — 00, which, when com- 

n— ^00 

bined with (jA4p . again gives liminf„^oo Ln(p) = 0, 
as required. Uniformity of convergence on any 
closed interval contained within (j)c, 1) follows as 
before, and, when taken with the earlier result, 
gives the uniformity assertion of the lemma. 

■ 

Proof. (Theorem IIII.ip We assume only that 
lim„_>oo in(p) exists. Then, by Lemma IIII.2I this limit 



is equal to zero for all p ^ Pc- (The existence of the limit 
would follow, for example, from the widely believed, but 
not rigorously proved, result that for p — Pc the distri- 
bution of the size of the open cluster size containing the 
origin has a power-law tail — see Chapter 9 of tlB].) Then 
the result of Lemma riII.2l can be reformulated in e-terms 
as follows: Ve>03iV(e,7)>0, such that 

£n{p) <eCn{Pc)yn> N{e,-f)\/pe [a,Pc-7MPc+l, l3), 

(A5) 

for any a and /?, such that 

A= [a,Pe-7]U[pc + 7,/3], Ac (0,1). 

The quantity Pn, being the maximum likelihood esti- 
mate for p, is the mode of £„(p): 

Pn := arg max -C„(p), 

P6(0,l) 

that is 

CniPn)>Cnip)'ipe{0,l). (A6) 

Consider the sequence of mle's {pn}'^=i- We now prove 
that this sequence converges to Pc- Suppose, conversely, 
this is not the case: 

3Ce (0,l)VAf > 037i> Af : \pn ~ Pc\ > C 

TakeM(C) = N{CX/2), then3n > M(C) : \pn-Pc\ > C 
and, thus, Pc ^ Pn- At the same time (when e = C) the 
following holds by \k5t\ : 

Cn{Pn) < C,Cn[Pc) < Cn{Pc): 

thereby contradicting (|A6p . Hence, Pn ^ Pc, n ^ 00. ■ 

Appendix B: Mixing properties of chains from the 
proposed algorithms 

Scenario 51 and Algorithm [H Figure [9] displays a 
trace plot of updates in p for inference made for the site 
configurations from Figure [T] (Algorithm [1] was applied) 
and presented in Figure [H The trace plot indicates that 
the mixing properties of the chain are rather satisfactory. 
The Matlab-based simulation used required 15 seconds 
on Intel(R) Core(TM)2 Duo CPU 2.26GHz to obtain a 
series of chain updates of the length 10^. 

Experiments with larger configurations suggest that 
the burn-in time of the constructed chain may vary con- 
siderably with the size of the site configuration and its 
potential connectivity properties as well as on the choice 
of the initial graph Go (in our experiments we worked 
with fully saturated graphs and with connected com- 
ponents derived from those by random sparsification of 
edges). However, in practice, the chain described by 
Algorithm [T] typically converges rapidly: both the con- 
nected component updates and the percolation param- 
eter updates can be efficiently realised, allowing one to 
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0.65 




2000 4000 6000 8000 10000 

step 

FIG. 9; Trace plot for McMC sampling resulted in the his- 
togram from Figure [2] for the cluster C in Figure [1] 



generate a suitable sample from the posterior distribution 
■nij) I C) within practically useful timescales. Updates of 
the connected component, that is deletion and insertion 
of edges while preserving the connectedness of the un- 
derlying graph, can be further improved using dynamic 
graph algorithms [13 ■ 

Scenario 52 and Algorithm[2l Figure [TU] shows his- 
tograms of samples from the distribution 7r(p| |C| = n) 
obtained using the proposed McMC for the scenario 52 
for the cluster size values p = 10, 35, 50, 70 and corre- 
sponding trace plots. 
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FIG. 10: Inference on the percolation parameter using McMC described in Algorithm (2] histograms of obtained samples and 
trace plots for (a,b) n = 10; (c,d) n — 35; (e,f) n = 50; (g,h) n = 70. 



